####################################
# Replication code for Study 1
## Paper: RACE, DIVERSITY, AND THE DEVELOPMENT OF POLITICAL ATTITUDES ON COLLEGE CAMPUSES
## Journal: JOURNAL OF POLITICS
## Authors: NATHAN KAR MING CHAN & TANIKA RAYCHAUDHURI 
####################################

##################
# Load packages #
## Note run using R: Version 4.3.1 or earlier
##################
rm(list=ls())

library(ggplot2)
library(gridExtra)
library(lme4)
library(lmerTest)
library(MatchIt)
library(stats)
library(texreg)


##################
# Load data #
##################

## set working directory -- replace with file path on your device
##setwd("")

## read in Study 1 data
data <- read.csv("study1data.csv", header = T) ## N = 401871 obs of 39 variables
ls(data)

##########################################
##########################################
## MAIN PAPER: STUDY 1
## Sample sizes listed in body of paper
##########################################
##########################################
## schools
schools <- unique(data$ACERECODE) ## 610 schools

## racial groups
table(data$RACEGROUP) ## 21,136 Asian; 24,029 Black, 19,242 Latino, 337,464 white

## years available 
summary(data$YEAR_TFS) ## 1987-2005: Freshman year
summary(data$YEAR) ## 1994-2006: Senior year


##########################################
##########################################
## MAIN PAPER: STUDY 1
## Code to create: Table 1; Figures 1-3
##########################################
##########################################

#######################################
## Main paper Table 1
## Changes in perceptions that racial discrimination is no longer 
## a major problem in the U.S. from freshman to senior year (CIRP panel survey)
## Note: fill in values based on code. Round to third decimal point.
#######################################

## All students (Table 1, Row 1)
tfs.all <- round(mean(data$discrimtfs, na.rm = T), 3) ## .259
css.all <- round(mean(data$discrimcss, na.rm = T),3) ## .229
diff.all <- round(mean(data$discrimcss, na.rm = T), 3) - round(mean(data$discrimtfs, na.rm = T),3) ##-0.030
t.test(data$discrimtfs, data$discrimcss, paired = F) ## p < 0.001

## African American students (Table 1, Row 2)
tfs.black <- round(mean(data$discrimtfs[data$black == 1], na.rm = T), 3) ## .119
css.black <- round(mean(data$discrimcss[data$black == 1], na.rm = T),3) ## .105
diff.black <- round(mean(data$discrimcss[data$black == 1], na.rm = T), 3) - round(mean(data$discrimtfs[data$black == 1], na.rm = T),3) ##-0.014
t.test(data$discrimtfs[data$black == 1], data$discrimcss[data$black == 1], paired = F) ## p < 0.001

## Asian American students 
tfs.asian <- round(mean(data$discrimtfs[data$asian == 1], na.rm = T), 3) ## .230
css.asian <- round(mean(data$discrimcss[data$asian == 1], na.rm = T),3) ## .201
diff.asian <- round(mean(data$discrimcss[data$asian == 1], na.rm = T), 3) - round(mean(data$discrimtfs[data$asian == 1], na.rm = T),3) ##-0.029
t.test(data$discrimtfs[data$asian == 1], data$discrimcss[data$asian == 1], paired = F) ## p < 0.001

## Latino students
tfs.latino <- round(mean(data$discrimtfs[data$latino == 1], na.rm = T), 3) ## .212
css.latino <- round(mean(data$discrimcss[data$latino == 1], na.rm = T),3) ## .179
diff.latino <- round(mean(data$discrimcss[data$latino == 1], na.rm = T), 3) - round(mean(data$discrimtfs[data$latino == 1], na.rm = T),3) ##-0.033
t.test(data$discrimtfs[data$latino == 1], data$discrimcss[data$latino == 1], paired = F) ## p < 0.001

## White students
tfs.white <- round(mean(data$discrimtfs[data$white == 1], na.rm = T), 3) ## .270
css.white <- round(mean(data$discrimcss[data$white == 1], na.rm = T),3) ## .243
diff.white <- round(mean(data$discrimcss[data$white == 1], na.rm = T), 3) - round(mean(data$discrimtfs[data$white == 1], na.rm = T),3) ##-0.027
t.test(data$discrimtfs[data$white == 1], data$discrimcss[data$white == 1], paired = F) ## p < 0.001


#######################################
## Association between ethnic studies and perceptions of discrimination
#######################################

#######################################
## Main paper: Figures 1 and 2
## Association between taking ethnic studies course and 
## changes in perceptions of racial discrimination (model with controls)
#######################################

########################################
## Unmatched models with controls
########################################

## All students

discrim.all.controls <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + ideotfs +
                      midwest + west + south + hbcu + college + asian + black + latino +
                      female + englishfl + protestant + jewish + othornorelig +
                      citizen + humantfs + sciencetfs + businesstfs + parincome +
                      + factor(YEAR) + (1| ACERECODE), data = data)
summary(discrim.all.controls)

## extract coefficients for plot
discrim.all.coefs <- data.frame(coef(summary(discrim.all.controls)))
discrim.all.coefs$p.z <- 2 * (1 - pnorm(abs(discrim.all.coefs$t.value)))
discrim.all.coefs

## White students
discrim.white.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                        midwest + west + south + hbcu + college +
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                        + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "White", ])
summary(discrim.white.controls)

## extract coefficients for plot
discrim.white.coefs <- data.frame(coef(summary(discrim.white.controls)))
discrim.white.coefs$p.z <- 2 * (1 - pnorm(abs(discrim.white.coefs$t.value)))
discrim.white.coefs

## Asian students
discrim.asian.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                        midwest + west + south + hbcu + college +
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                        + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Asian", ])
summary(discrim.asian.controls)

## extract coefficients for plot
discrim.asian.coefs <- data.frame(coef(summary(discrim.asian.controls)))
discrim.asian.coefs$p.z <- 2 * (1 - pnorm(abs(discrim.asian.coefs$t.value)))
discrim.asian.coefs

## Latino students
discrim.latino.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                         midwest + west + south + hbcu + college +
                         female + englishfl + protestant + jewish + othornorelig +
                         citizen + humantfs + sciencetfs + businesstfs + parincome +
                         + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Hispanic", ])
summary(discrim.latino.controls)

## extract coefficients for plot
discrim.latino.coefs <- data.frame(coef(summary(discrim.latino.controls)))
discrim.latino.coefs$p.z <- 2 * (1 - pnorm(abs(discrim.latino.coefs$t.value)))
discrim.latino.coefs

## Black students
discrim.black.controls<- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                        midwest + west + south + hbcu + college +
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                        + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Black", ])
summary(discrim.black.controls)

## extract coefficients for plot
discrim.black.coefs <- data.frame(coef(summary(discrim.black.controls)))
discrim.black.coefs$p.z <- 2 * (1 - pnorm(abs(discrim.black.coefs$t.value)))
discrim.black.coefs

#######################################
## matched models w. controls
## full sample
#######################################

## create version of dataset without missing values
data.cem <- data[!is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
matching1 <- matchit(ethnicstudies ~ ideotfs + female + citizen + white + englishfl + south + hbcu + college, data = data.cem, 
                     method = 'cem', estimand = 'ATE')
summary(matching1, un=FALSE)

# Extract matched data
matched_df1 <- match.data(matching1)
#View(matched_df1)

discrim.matched <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + factor(YEAR) + (1| ACERECODE), data = matched_df1)
summary(discrim.matched)


discrim.matched.control <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + ideotfs +
                                  midwest + west + south + hbcu + college + asian + black + latino +
                                  female + englishfl + protestant + jewish + othornorelig +
                                  citizen + humantfs + sciencetfs + businesstfs + parincome +
                                  + factor(YEAR) + (1| ACERECODE), data = matched_df1)
summary(discrim.matched.control)

## extract coefficients for plot
discrim.matched.coefs <- data.frame(coef(summary(discrim.matched.control)))
discrim.matched.coefs


#######################################
## matched models - White
#######################################

## create version of dataset without missing values
data.white.cem <- data[data$white == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
matching3 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.white.cem, 
                     method = 'cem', estimand = 'ATE')
summary(matching3, un=FALSE)

# Extract matched data
matched_df3 <- match.data(matching3)
##View(matched_df3) 

## White
discrim.matched.control.white <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                                        midwest + west + south + hbcu + college +
                                        female + englishfl + protestant + jewish + othornorelig +
                                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                                        + factor(YEAR) + (1| ACERECODE), data = matched_df3)
summary(discrim.matched.control.white)

## extract coefficients for plot
discrim.matched.white.coefs <- data.frame(coef(summary(discrim.matched.control.white)))
discrim.matched.white.coefs


#######################################
## matched models - Asian
#######################################
## Asian
## create version of dataset without missing values
data.asian.cem <- data[data$asian == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
matching2 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.asian.cem, 
                     method = 'cem', estimand = 'ATE')
summary(matching2, un=FALSE)

# Extract matched data
matched_df2 <- match.data(matching2)
#View(matched_df2)

discrim.matched.control.asian <- lmer(discrimchg ~  ethnicstudies + racialorg +  racialws + ideotfs +
                                        midwest + west + south + hbcu + college +
                                        female + englishfl + protestant + jewish + othornorelig +
                                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                                        + factor(YEAR) + (1| ACERECODE), data = matched_df2)
summary(discrim.matched.control.asian)

## extract coefficients for plot
discrim.matched.asian.coefs <- data.frame(coef(summary(discrim.matched.control.asian)))
discrim.matched.asian.coefs

#######################################
## matched models - Latino
#######################################

## create version of dataset without missing values
data.latino.cem <- data[data$latino == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
matching5 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.latino.cem, 
                     method = 'cem', estimand = 'ATE')
summary(matching5, un=FALSE)

# Extract matched data
matched_df5 <- match.data(matching5)
#View(matched_df5)

discrim.matched.control.latino <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                                         midwest + west + south + hbcu + college +
                                         female + englishfl + protestant + jewish + othornorelig +
                                         citizen + humantfs + sciencetfs + businesstfs + parincome +
                                         + factor(YEAR) + (1| ACERECODE), data = matched_df5)
summary(discrim.matched.control.latino)

## extract coefficients for plot
discrim.matched.latino.coefs <- data.frame(coef(summary(discrim.matched.control.latino)))
discrim.matched.latino.coefs


#######################################
## matched models - Black
#######################################
## Black
## create version of dataset without missing values
data.black.cem <- data[data$black == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
matching4 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.black.cem, 
                     method = 'cem', estimand = 'ATE')
summary(matching4, un=FALSE)

# Extract matched data
matched_df4 <- match.data(matching4)
#View(matched_df4)

discrim.matched.control.black <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                                        midwest + west + south + hbcu + college +
                                        female + englishfl + protestant + jewish + othornorelig +
                                        citizen + humantfs + sciencetfs + businesstfs + parincome +
                                        + factor(YEAR) + (1| ACERECODE), data = matched_df4)
summary(discrim.matched.control.black)

## extract coefficients for plot
discrim.matched.black.coefs <- data.frame(coef(summary(discrim.matched.control.black)))
discrim.matched.black.coefs


#######################################
## Plot - models with controls
#######################################


data2 = data.frame (Model = c("Standard HLM", 
                              "Matched HLM"),
                    effect = c(discrim.all.coefs$Estimate[2], discrim.matched.coefs$Estimate[2],
                               discrim.white.coefs$Estimate[2], discrim.matched.white.coefs$Estimate[2],
                               discrim.asian.coefs$Estimate[2], discrim.matched.asian.coefs$Estimate[2],
                               discrim.black.coefs$Estimate[2], discrim.matched.black.coefs$Estimate[2],
                               discrim.latino.coefs$Estimate[2], discrim.matched.latino.coefs$Estimate[2]),
                    stderror = c(discrim.all.coefs$Std..Error[2], discrim.matched.coefs$Std..Error[2],
                                 discrim.white.coefs$Std..Error[2], discrim.matched.white.coefs$Std..Error[2],
                                 discrim.asian.coefs$Std..Error[2], discrim.matched.asian.coefs$Std..Error[2],
                                 discrim.black.coefs$Std..Error[2], discrim.matched.black.coefs$Std..Error[2],
                                 discrim.latino.coefs$Std..Error[2], discrim.matched.latino.coefs$Std..Error[2]),
                    Race = c("All respondents", "All respondents",
                             "White respondents", "White respondents",
                             "Asian respondents", "Asian respondents",
                             "Black respondents", "Black respondents",
                             "Latino/a respondents", "Latino/a respondents")) 

data2$Model <- factor(data2$Model, levels = c("Standard HLM", 
                                              "Matched HLM"))

data2 


p1.all <- ggplot(data2[1:2,], aes(x=Model, y=effect, linetype = Model)) +
  facet_grid(.~Race) +
  geom_point(stat="identity") +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3) +
  ylab("Change in perception of racial discrimination") +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), legend.position = "none") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.025, 0.005))
p1.all

p2.white <- ggplot(data2[3:4,], aes(x=Model, y=effect, linetype = Model)) +
  facet_grid(.~Race) +
  geom_point(stat="identity") +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3) +
  #ylab("Change in perception that racial discrimination is no longer a problem") +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), axis.title.y=element_blank(), legend.position = "none") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.045, 0.015))
p2.white

p3.asian <- ggplot(data2[5:6,], aes(x=Model, y=effect, linetype = Model)) +
  facet_grid(.~Race) +
  geom_point(stat="identity") +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3) +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), axis.title.y=element_blank(),legend.position = "none") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.045, 0.015))
p3.asian

p4.black <- ggplot(data2[7:8,], aes(x=Model, y=effect, linetype = Model)) +
  facet_grid(.~Race) +
  geom_point(stat="identity") +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3) +
  #ylab("Change in perception that racial discrimination is no longer a problem") +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), axis.title.y=element_blank(), legend.position = "none") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.045, 0.015))
p4.black

p5.latino <- ggplot(data2[9:10,], aes(x=Model, y=effect, linetype = Model)) +
  facet_grid(.~Race) +
  geom_point(stat="identity") +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3) +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), axis.title.y=element_blank(), legend.position = "none") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.045, 0.015))
p5.latino

#######################################
## Figure 1: Association between taking ethnic studies 
## course and changes in perceptions of racial discrimination (model with controls)
#######################################

pdf("figure1.pdf", width = 5, height = 6)
grid.arrange(p1.all)
dev.off()

#######################################
## Figure 2: Association between taking an ethnic studies 
## course and changes in perceptions of racial discrimination
## (disaggregated across racial groups with controls)
#######################################

pdf("figure2.pdf", width = 9, height = 6)
grid.arrange(p2.white, p3.asian, p4.black, p5.latino, nrow = 2,
             left = "Change in perceptions of racial discrimination")
dev.off()

#######################################
## Selection test: Perceptions of discrimination 
#######################################

## Code for Main Paper text: "In these data, students who enrolled in ethnic studies courses entered college 
## with lower perceptions that discrimination is no longer a problem than those who did not enroll in such courses." 
## mean levels for those who took and did not take ethnic studies"
tapply(data$discrimtfs, data$ethnicstudies, mean, na.rm = T)
t.test(data$discrimtfs[data$ethnicstudies == 1], data$discrimtfs[data$ethnicstudies == 0], paired = F)

## create terciles and subset data
table(data$discrimtfs)
quantile(data$discrimtfs, c(0.33, 0.66, 1), na.rm = T)

data$lowdiscrimtfs <- ifelse(data$discrimtfs == 0, 1, 0)
table(data$lowdiscrimtfs) ## 88,388
data$middiscrimtfs <- ifelse(data$discrimtfs > 0 & data$discrimtfs < 0.34, 1, 0)
table(data$middiscrimtfs)  ## 97,446 
data$highdiscrimtfs <- ifelse(data$discrimtfs > 0.33, 1, 0)
table(data$highdiscrimtfs)  ## 132,002 

data.lowdisc <- data[data$lowdiscrimtfs == 1 & !is.na(data$lowdiscrimtfs == 1), ]  ## 88,388 obs
data.middisc <- data[data$middiscrimtfs == 1 & !is.na(data$middiscrimtfs == 1), ] ## 97,446 obs
data.highdisc <- data[data$highdiscrimtfs == 1 & !is.na(data$highdiscrimtfs == 1), ] ## 132,002 obs

## basic models
discrim.all.low <- lmer(discrimchg ~ ethnicstudies + racialorg + 
                           racialws + factor(YEAR) + (1| ACERECODE), data = data.lowdisc)
summary(discrim.all.low)

## extract coefficients for plot
discrim.low.coefs <- data.frame(coef(summary(discrim.all.low)))
discrim.low.coefs


discrim.all.mid <- lmer(discrimchg ~ ethnicstudies + racialorg + 
                           racialws + factor(YEAR) + (1| ACERECODE), data = data.middisc)
summary(discrim.all.mid)

## extract coefficients for plot
discrim.mid.coefs <- data.frame(coef(summary(discrim.all.mid)))
discrim.mid.coefs

discrim.all.high <- lmer(discrimchg ~ ethnicstudies +  racialorg + 
                           racialws + factor(YEAR) + (1| ACERECODE), data = data.highdisc)
summary(discrim.all.high)

## extract coefficients for plot
discrim.high.coefs <- data.frame(coef(summary(discrim.all.high)))
discrim.high.coefs

## models with controls
discrim.all.low.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                               midwest + west + south + hbcu + college + asian + black + latino +
                               female + englishfl + protestant + jewish + othornorelig +
                               citizen + humantfs + sciencetfs + businesstfs + parincome +
                               + factor(YEAR) + (1| ACERECODE), data = data.lowdisc)
summary(discrim.all.low.cont)

## extract coefficients for plot
discrim.low.cont.coefs <- data.frame(coef(summary(discrim.all.low.cont)))
discrim.low.cont.coefs

discrim.all.mid.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                               midwest + west + south + hbcu + college + asian + black + latino +
                               female + englishfl + protestant + jewish + othornorelig +
                               citizen + humantfs + sciencetfs + businesstfs + parincome +
                               + factor(YEAR) + (1| ACERECODE), data = data.middisc)
summary(discrim.all.mid.cont)

## extract coefficients for plot
discrim.mid.cont.coefs <- data.frame(coef(summary(discrim.all.mid.cont)))
discrim.mid.cont.coefs

discrim.all.high.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                                midwest + west + south + hbcu + college + asian + black + latino +
                                female + englishfl + protestant + jewish + othornorelig +
                                citizen + humantfs + sciencetfs + businesstfs + parincome +
                                + factor(YEAR) + (1| ACERECODE), data = data.highdisc)
summary(discrim.all.high.cont)

## extract coefficients for plot
discrim.high.cont.coefs <- data.frame(coef(summary(discrim.all.high.cont)))
discrim.high.cont.coefs



#######################################
## Plot - pre-college discrim perception plots
#######################################


data3 = data.frame (Model = c("Basic HLM", 
                              "Controlled HLM"),
                    effect = c(discrim.low.coefs$Estimate[2], discrim.low.cont.coefs$Estimate[2],
                               discrim.mid.coefs$Estimate[2], discrim.mid.cont.coefs$Estimate[2],
                               discrim.high.coefs$Estimate[2], discrim.high.cont.coefs$Estimate[2]),
                    stderror = c(discrim.low.coefs$Std..Error[2], discrim.low.cont.coefs$Std..Error[2],
                                 discrim.mid.coefs$Std..Error[2], discrim.mid.cont.coefs$Std..Error[2],
                                 discrim.high.coefs$Std..Error[2], discrim.high.cont.coefs$Std..Error[2]),
                    Discrimination = c("Low",  "Low", 
                                       "Medium", "Medium",  
                                       "High", "High")) 

data3$Discrimination <- factor(data3$Discrimination, levels = c("Low", "Medium", "High"))


data3

p1.discrim <- ggplot(data3, aes(x=Model, y=effect, linetype = Discrimination)) +  geom_point(stat="identity", position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(ymin = (effect) - 1.96*(stderror), 
                    ymax = (effect) + 1.96*(stderror)), width = .3, position = position_dodge(width = 0.9)) +
  ylab("Change in perceptions of racial discrmination") +
  #guides(fill=FALSE) +
  theme(text = element_text(size=12), axis.title.x=element_blank(), legend.position = "bottom") +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(limits=c(-0.042, 0.001))+ 
  labs(color = "Pre-college perceptions of discrimination")
p1.discrim

#######################################
## Figure 3: Association between taking an ethnic studies 
## course and changes in perceptions of racial discrimination
## (conditioning on pre-college perceptions of discrimination with controls)
#######################################

pdf("figure3.pdf", width = 7.5, height = 6)
grid.arrange(p1.discrim, nrow = 1)
dev.off()

##########################################
##########################################
## APPENDIX: STUDY 1 
## Code to create Appendix Tables A1-A10
##########################################
##########################################

##########################################
## Appendix Table A1: (Regression output)
## Association between enrolling in ethnic studies course and 
## changes in perceptions of racial discrimination (Unmatched sample, no controls)
##########################################

#######################################
## Unmatched models
## Full sample and all racial groups
#######################################

## All students
discrim.all <- lmer(discrimchg ~ ethnicstudies + racialorg + 
                      racialws + factor(YEAR) + (1| ACERECODE), data = data)
summary(discrim.all)


## White students
discrim.white <- lmer(discrimchg ~  ethnicstudies +racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "White", ])
summary(discrim.white)


## Asian students
discrim.asian <- lmer(discrimchg ~  ethnicstudies + racialorg +  racialws + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Asian", ])
summary(discrim.asian)


## Latino students
discrim.latino <- lmer(discrimchg ~   ethnicstudies + racialorg +  racialws + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Hispanic", ])
summary(discrim.latino)

## Black students
discrim.black <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Black", ])
summary(discrim.black)

## Print Table A1:
htmlreg(list(discrim.all, discrim.white, discrim.asian, discrim.latino, discrim.black),
        custom.coef.names= c("Intercept","Enrolled in ethnic studies course", "Participated in an ethnic or racial student org",
                             "Attended racial awareness workshop", "Year: 1995", "Year: 1996", "Year: 1997", "Year: 1998", "Year: 1999",
                             "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Full sample", "White", "Asian", "Latino", "Black"),
        digits = 3,
        file = "TableA1.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients. 

##########################################
## Appendix Table A2: (Regression output)
## Association between enrolling in ethnic studies course and 
## changes in perceptions of racial discrimination (Matched sample, no controls)
##########################################

#######################################
## matched models- Full sample (no controls)
## Note: uses matched dataframe created earlier (code included if needed)
#######################################
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
## create version of dataset without missing values
#data.cem <- data[!is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
#matching1 <- matchit(ethnicstudies ~ ideotfs + female + citizen + white + englishfl + south + hbcu + college, data = data.cem, 
#                     method = 'cem', estimand = 'ATE')
#summary(matching1, un=FALSE)

# Extract matched data
#matched_df1 <- match.data(matching1)
#View(matched_df1)

discrim.matched <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + factor(YEAR) + (1| ACERECODE), data = matched_df1)
summary(discrim.matched)

#######################################
## matched models - White (no controls)
## Note: uses matched dataframe created earlier (code included if needed)
#######################################
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
## create version of dataset without missing values
#data.white.cem <- data[data$white == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
#matching3 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.white.cem, 
#                     method = 'cem', estimand = 'ATE')
#summary(matching3, un=FALSE)

# Extract matched data
#matched_df3 <- match.data(matching3)
##View(matched_df3)

discrim.matched.white <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = matched_df3)
summary(discrim.matched.white)

#######################################
## matched models - Asian
## Note: uses matched dataframe created earlier (code included if needed)
#######################################
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
## create version of dataset without missing values
#data.asian.cem <- data[data$asian == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
#matching2 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.asian.cem, 
#                     method = 'cem', estimand = 'ATE')
#summary(matching2, un=FALSE)

# Extract matched data
#matched_df2 <- match.data(matching2)
#View(matched_df2)

discrim.matched.asian <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = matched_df2)
summary(discrim.matched.asian)


#######################################
## matched models - Latino (no controls)
## Note: uses matched dataframe created earlier (code included if needed)
#######################################
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
## create version of dataset without missing values
#data.latino.cem <- data[data$latino == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
#matching5 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.latino.cem, 
#                     method = 'cem', estimand = 'ATE')
#summary(matching5, un=FALSE)

# Extract matched data
#matched_df5 <- match.data(matching5)
#View(matched_df5)

discrim.matched.latino <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = matched_df5)
summary(discrim.matched.latino)


#######################################
## matched models - Black (no controls)
## Note: uses matched dataframe created earlier (code included if needed)
#######################################
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
## create version of dataset without missing values
#data.black.cem <- data[data$black == 1 & !is.na(data$ideotfs), ]

# Perform CEM --------------------
## match on ideology, gender, race citizenship, englishfl, south, hbcu, college
#matching4 <- matchit(ethnicstudies ~ ideotfs + female + citizen + englishfl + south + hbcu + college, data = data.black.cem, 
#                     method = 'cem', estimand = 'ATE')
#summary(matching4, un=FALSE)

# Extract matched data
#matched_df4 <- match.data(matching4)
#View(matched_df4)

discrim.matched.black <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + factor(YEAR) + (1| ACERECODE), data = matched_df4)
summary(discrim.matched.black)

## Print Table A2:
htmlreg(list(discrim.matched, discrim.matched.white, discrim.matched.asian, discrim.matched.latino, discrim.matched.black),
        custom.coef.names= c("Intercept","Enrolled in ethnic studies course", "Participated in an ethnic or racial student org",
                             "Attended racial awareness workshop", "Year: 1995", "Year: 1996", "Year: 1997", "Year: 1998", "Year: 1999",
                             "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Full sample", "White", "Asian", "Latino", "Black"),
        digits = 3,
        file = "TableA2.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients. 

##########################################
## Appendix Table A3: (Regression output)
## Association between enrolling in ethnic studies course and 
## changes in perceptions of racial discrimination (Unmatched sample, with controls)
## Prints regression models included above. (Code included if needed)
##########################################

## All students
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.all.controls <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + ideotfs +
#                               midwest + west + south + hbcu + college + asian + black + latino +
#                               female + englishfl + protestant + jewish + othornorelig +
#                               citizen + humantfs + sciencetfs + businesstfs + parincome +
#                               + factor(YEAR) + (1| ACERECODE), data = data)
summary(discrim.all.controls)

## White 
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.white.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
#                                 midwest + west + south + hbcu + college +
#                                 female + englishfl + protestant + jewish + othornorelig +
#                                 citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                 + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "White", ])
summary(discrim.white.controls)


## Asian
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.asian.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
#                                 midwest + west + south + hbcu + college +
#                                 female + englishfl + protestant + jewish + othornorelig +
#                                 citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                 + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Asian", ])
summary(discrim.asian.controls)


## Latino 
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.latino.controls <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
#                                  midwest + west + south + hbcu + college +
#                                  female + englishfl + protestant + jewish + othornorelig +
#                                  citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                  + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Hispanic", ])
summary(discrim.latino.controls)


## Black
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.black.controls<- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                                midwest + west + south + hbcu + college +
#                                female + englishfl + protestant + jewish + othornorelig +
#                                citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                + factor(YEAR) + (1| ACERECODE), data = data[data$RACEGROUP == "Black", ])
summary(discrim.black.controls)

## Print Table A3:
htmlreg(list(discrim.all.controls, discrim.white.controls, discrim.asian.controls, discrim.latino.controls, discrim.black.controls),
        custom.coef.names = c("Intercept", "Enrolled in ethnic studies course", "Participated in an ethnic or racial student org",
                              "Attended racial/cultural awareness workshop", 
                              "Ideology", "Midwest", "West", "South", "HBCU", "College", "Asian", "Black", "Latino",
                              "Female", "English first language", "Protestant", "Jewish", "Other or no religion",
                              "Citizen", "Humanities", "Science", "Business", "Parental income", "Year: 1995", "Year: 1996", "Year: 1997",
                              "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Full sample", "White", "Asian", "Latino", "Black"),
        digits = 3,
        file = "TableA3.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients. 

##########################################
## Appendix Table A4: (Regression output)
## Association between enrolling in ethnic studies course and 
## changes in perceptions of racial discrimination (Matched sample, with controls)
## Prints regression models included above. (Code included if needed)
##########################################
## All students
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
#discrim.matched.control <- lmer(discrimchg ~  ethnicstudies + racialorg + racialws + ideotfs +
#                                  midwest + west + south + hbcu + college + asian + black + latino +
#                                  female + englishfl + protestant + jewish + othornorelig +
#                                  citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                  + factor(YEAR) + (1| ACERECODE), data = matched_df1)
summary(discrim.matched.control)

## White
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
#discrim.matched.control.white <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                                        midwest + west + south + hbcu + college +
#                                        female + englishfl + protestant + jewish + othornorelig +
#                                        citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                        + factor(YEAR) + (1| ACERECODE), data = matched_df3)
summary(discrim.matched.control.white)

## Asian
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
#discrim.matched.control.asian <- lmer(discrimchg ~  ethnicstudies + racialorg +  racialws + ideotfs +
#                                        midwest + west + south + hbcu + college +
#                                       female + englishfl + protestant + jewish + othornorelig +
#                                        citizen + humantfs + sciencetfs + businesstfs + parincome +                                        
#                                       + factor(YEAR) + (1| ACERECODE), data = matched_df2)
summary(discrim.matched.control.asian)

## Latino
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
#discrim.matched.control.latino <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
#                                         midwest + west + south + hbcu + college +
#                                         female + englishfl + protestant + jewish + othornorelig +
#                                         citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                         + factor(YEAR) + (1| ACERECODE), data = matched_df5)

summary(discrim.matched.control.latino)

## Black
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.
#discrim.matched.control.black <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                                        midwest + west + south + hbcu + college +
#                                        female + englishfl + protestant + jewish + othornorelig +
#                                        citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                        + factor(YEAR) + (1| ACERECODE), data = matched_df4)
summary(discrim.matched.control.black)

## Print Table A4:
htmlreg(list(discrim.matched.control, discrim.matched.control.white, discrim.matched.control.asian, 
             discrim.matched.control.latino, discrim.matched.control.black),
        custom.coef.names = c("Intercept", "Enrolled in ethnic studies course", "Participated in an ethnic or racial student org",
                              "Attended racial/cultural awareness workshop", 
                              "Ideology", "Midwest", "West", "South", "HBCU", "College", 
                              "Asian", "Black", "Latino",
                              "Female", "English first language", "Protestant", "Jewish", "Other or no religion",
                              "Citizen", "Humanities", "Science", "Business", "Parental income", "Year: 1995", "Year: 1996", "Year: 1997",
                              "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Full sample", "White", "Asian", "Latino", "Black"),
        digits = 3,
        file = "TableA4.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients.

##########################################
## Appendix Table A5: (Regression output)
## Association between enrolling in ethnic studies course and changes in 
## perceptions of racial discrimination (disaggregated by pre-college perceptions of discrimination, with controls)
## Prints regression models included above. (Code included if needed)
##########################################

## Models without controls 
## Note: Uncomment the following lines if running this section of code separately from Main paper figures.

#discrim.all.low <- lmer(discrimchg ~  ethnicstudies + racialorg + 
#                          racialws + factor(YEAR) + (1| ACERECODE), data = data.lowdisc)
summary(discrim.all.low)


#discrim.all.mid <- lmer(discrimchg ~ ethnicstudies + racialorg + 
#                          racialws + factor(YEAR) + (1| ACERECODE), data = data.middisc)
summary(discrim.all.mid)


#discrim.all.high <- lmer(discrimchg ~ ethnicstudies + racialorg + 
#                           racialws + factor(YEAR) + (1| ACERECODE), data = data.highdisc)
summary(discrim.all.high)

## models with controls
#discrim.all.low.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                               midwest + west + south + hbcu + college + asian + black + latino +
#                               female + englishfl + protestant + jewish + othornorelig +
#                               citizen + humantfs + sciencetfs + businesstfs + parincome +
#                               + factor(YEAR) + (1| ACERECODE), data = data.lowdisc)
summary(discrim.all.low.cont)

#discrim.all.mid.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                               midwest + west + south + hbcu + college + asian + black + latino +
#                               female + englishfl + protestant + jewish + othornorelig +
#                               citizen + humantfs + sciencetfs + businesstfs + parincome +
#                               + factor(YEAR) + (1| ACERECODE), data = data.middisc)
summary(discrim.all.mid.cont)


#discrim.all.high.cont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
#                                midwest + west + south + hbcu + college + asian + black + latino +
#                                female + englishfl + protestant + jewish + othornorelig +
#                                citizen + humantfs + sciencetfs + businesstfs + parincome +
#                                + factor(YEAR) + (1| ACERECODE), data = data.highdisc)
summary(discrim.all.high.cont)

## Print Table A5:
htmlreg(list(discrim.all.high, discrim.all.high.cont,
             discrim.all.mid, discrim.all.mid.cont, 
             discrim.all.low, discrim.all.low.cont),
        custom.coef.names = c("Intercept", "Enrolled in ethnic studies course",  "Participated in an ethnic or racial student org",
                              "Attended racial/cultural awareness workshop", "Year: 1995", "Year: 1996", "Year: 1997",
                              "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006",
                              "Ideology", "Midwest", "West", "South", "HBCU", "College", "Asian", "Black", "Latino",
                              "Female", "English first language", "Protestant", "Jewish", "Other or no religion",
                              "Citizen", "Humanities", "Science", "Business", "Parental income"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("High perception", "High perception", 
                               "Medium perception", "Medium perception",  
                               "Low perception", "Low perception"),
        digits = 3,
        file = "TableA5.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients.

##########################################
## Appendix Table A6: (Regression output)
## Association between enrolling in ethnic studies course and 
## changes in perceptions of racial discrimination 
## (Unmatched sample, with cohort racial demographic controls)
##########################################

discrim.all.newcont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                      midwest + west + south + hbcu + college + asian + black + latino +
                      female + englishfl + protestant + jewish + othornorelig +
                      citizen + humantfs + sciencetfs + businesstfs + parincome + 
                      blackper + asianper + latinoper + othraceper
                    + factor(YEAR) + (1| ACERECODE), data = data)
summary(discrim.all.newcont)


discrim.white.newcont <- lmer(discrimchg ~ ethnicstudies + racialorg +  racialws + ideotfs +
                        midwest + west + south + hbcu + college + 
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome + 
                        blackper + asianper + latinoper + othraceper
                      + factor(YEAR) + (1| ACERECODE), data = data[data$white == 1, ])
summary(discrim.white.newcont)

discrim.asian.newcont <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                        midwest + west + south + hbcu + college + 
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome + 
                        blackper + asianper + latinoper + othraceper
                      + factor(YEAR) + (1| ACERECODE), data = data[data$asian == 1, ])
summary(discrim.asian.newcont)

discrim.latino.newcont <- lmer(discrimchg ~ ethnicstudies + racialorg + racialws + ideotfs +
                         midwest + west + south + hbcu + college + asian + black + latino +
                         female + englishfl + protestant + jewish + othornorelig +
                         citizen + humantfs + sciencetfs + businesstfs + parincome + 
                         blackper + asianper + latinoper + othraceper
                       + factor(YEAR) + (1| ACERECODE), data = data[data$latino == 1, ])
summary(discrim.latino.newcont)


discrim.black.newcont <- lmer(discrimchg ~  ethnicstudies +  racialorg +racialws + ideotfs +
                        midwest + west + south + hbcu + college + 
                        female + englishfl + protestant + jewish + othornorelig +
                        citizen + humantfs + sciencetfs + businesstfs + parincome + 
                        blackper + asianper + latinoper + othraceper
                      + factor(YEAR) + (1| ACERECODE), data = data[data$black == 1, ])
summary(discrim.black.newcont)

## Print Table A6:
htmlreg(list(discrim.all.newcont, discrim.white.newcont, discrim.asian.newcont, discrim.latino.newcont, discrim.black.newcont),
        custom.coef.names = c("Intercept", "Enrolled in ethnic studies course", "Participated in an ethnic or racial student org",
                              "Attended racial/cutural awareness workshop", 
                              "Ideology", "Midwest", "West", "South", "HBCU", "College", "Asian", "Black", "Latino",
                              "Female", "English first language", "Protestant", "Jewish", "Other or no religion",
                              "Citizen", "Humanities", "Science", "Business", "Parental income", 
                              "Perc. Black", "Perc. Asian", "Perc. Latino", "Perc. Other race",
                              "Year: 1995", "Year: 1996", "Year: 1997",
                              "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Full sample", "White", "Asian", "Latino", "Black"),
        digits = 3,
        file = "TableA6.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients.

#######################################
## Appendix Table A7: 
## Institutional exposure to diversity correlation matrix
#######################################
deidata <- data[c("ethnicstudies", "racialws", "racialorg")]
cor_data = round(cor(deidata), 2)
cor_data ## Prints values in Table A7


#######################################
## Appendix Table A8: 
## Association between taking an ethnic studies course and 
## other measures of institutional exposure to diversity
#######################################

deiactivity.1 <- lmer(racialws ~ 
                        ethnicstudies + factor(YEAR) + (1| ACERECODE), data = data)
summary(deiactivity.1)

deiactivity.2 <- lmer(racialorg ~ 
                        ethnicstudies + factor(YEAR) + (1| ACERECODE), data = data)
summary(deiactivity.2)

## regression output
htmlreg(list(deiactivity.1, deiactivity.2),
        custom.coef.names = c ("Intercept", "Enrolled in ethnic studies course", "Year: 1995", "Year: 1996", "Year: 1997",
                               "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        star.cutoffs = c(.05, .01, .001),
        custom.model.names = c("Attended racial/cultural awareness workshop", "Participated in ethnic or racial student org"),
        digits = 3,
        file = "TableA8.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients.

#######################################
## Appendix Table A9: 
#######################################

ethnicstudies.individual <- lmer(ethnicstudies ~ black + asian + latino + female + ideotfs +
                                   citizen + parincome + humantfs + sciencetfs + businesstfs
                                 + factor(YEAR) + (1| ACERECODE), data = data)
summary(ethnicstudies.individual)

ethnicstudies.school <- lmer(ethnicstudies ~ midwest + south + west + hbcu + college 
                             + factor(YEAR) + (1| ACERECODE), data = data)
summary(ethnicstudies.school)

## regression output
htmlreg(list(ethnicstudies.individual, ethnicstudies.school),
        custom.coef.names = c ("Intercept", "African American", "Asian", "Latino", "Female", "Ideology", "Citizen",
                              "Parental income", "Humanities", "Science", "Buisness",
                              "Year: 1995", "Year: 1996", "Year: 1997", "Year: 1998", "Year: 1999", "Year: 2000", "Year: 2001", 
                              "Year: 2002", "Year: 2003", "Year: 2004", "Year: 2005", "Year: 2006",
                              "Midwest", "South", "West", "HBCU", "College"),
        omit.coef = "YEAR",
        star.char = c("*", "**", "***"),
        custom.model.names = c("Individual-level factors", "School-level factors"),
        star.cutoffs = c(.05, .01, .001),
        digits = 3,
        file = "TableA9.html")

## Notes: Final table in paper does not report AIC, BIC, Log Likelihood or Variance. 
## ACERECODE refers to the number of schools.
## Suppressed year fixed effect coefficients.